#include "AnaDef.h"
#include "XGLUtils.h"
#include "style.h"

#define EPSILON 1e-12

TH2D *fmcmom=0x0;
TH2D *fmcpt=0x0;
TH2D *fmcz=0x0;
TH2D *fmcxi=0x0;

void getMCETAall()
{
  const Int_t  grid[]={50,60,70,80, 53, 63,73, 83};
  TH2D * hmcs[]={fmcmom, fmcpt, fmcz, fmcxi, fmcmom, fmcpt, fmcz, fmcxi};
  const Int_t nid = sizeof(grid)/sizeof(Int_t);
  
  TString pwd = gSystem->pwd();
  pwd=pwd(pwd.Last('/')+1,10000);
  printf("pwd %s\n", pwd.Data());

  for(Int_t iid=0; iid<nid; iid++){
    TList *ll = new TList;
    ll->SetOwner(kTRUE);

    TH1D * hmcnor = 0x0;
    if(grid[iid]%10==0){
      hmcnor = hmcs[iid]->ProjectionX(Form("%sall", hmcs[iid]->GetName()), hmcs[iid]->GetYaxis()->FindBin(0.+EPSILON), hmcs[iid]->GetYaxis()->FindBin(3.+EPSILON));
    }
    else{
      hmcnor = hmcs[iid]->ProjectionX(Form("%spi",  hmcs[iid]->GetName()), hmcs[iid]->GetYaxis()->FindBin(1.+EPSILON), hmcs[iid]->GetYaxis()->FindBin(1.+EPSILON));
    }

    for(Int_t itype=0; itype<4; itype++){
      const Int_t tmpy = hmcs[iid]->GetYaxis()->FindBin(itype+EPSILON);
      TH1D * htype = hmcs[iid]->ProjectionX(Form("%sp%d0", hmcs[iid]->GetName(), itype), tmpy, tmpy);

      //htype->Divide(htype, hmcnor, 1, 1, "B");
      for(Int_t ibin=1; ibin<=htype->GetNbinsX(); ibin++){
        const Double_t ydown = hmcnor->GetBinContent(ibin);
        const Double_t yup = htype->GetBinContent(ibin);

        if(ydown<EPSILON || yup<EPSILON){
          htype->SetBinContent(ibin,0); htype->SetBinError(ibin,0);
          continue;
        }

        const Double_t rat = yup/ydown;
        Double_t err = 0;
        if(grid[iid]%10==0){
          err = TMath::Sqrt(rat*(1-rat)/ydown);
        }
        else{
          err = rat * TMath::Sqrt(1/yup+1/ydown); //correct since yup and ydown are both real counts
        }
        htype->SetBinContent(ibin, rat); htype->SetBinError(ibin,err);
      }

      TGraphAsymmErrors * gtype = XGLUtils::HistToGraph(Form("grstat%d", itype), htype); 
      ll->Add(gtype);
      if(grid[iid]/10<=6){
        style::ToNaturalScale(gtype);
      }
      TGraphAsymmErrors * gsys = (TGraphAsymmErrors*) gtype->Clone(Form("grsys%d", itype));
      for(Int_t ip=0; ip<gsys->GetN(); ip++){
        gsys->SetPointEYlow(ip,0);
        gsys->SetPointEYhigh(ip,0);
      }
      ll->Add(gsys);
      
      delete htype;
    }

    TFile *fout=new TFile(Form("%s_directMC_id%d.root", pwd.Data(), grid[iid]),"recreate");
    ll->Write();
    fout->Save();
    fout->Close();

    delete fout;
    delete ll;
    delete hmcnor;
  }
}

int main(int argc, char * argv[])
{
  TH1::AddDirectory(kFALSE);

  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=2){
    printf("drawpar argc!=2\n"); 
    return 1;
  }

  Int_t nrbx = 1;

  const TString pwd = gSystem->pwd();
  if(pwd.Contains("_rb3")){
    nrbx = 3;
  }
  else if(pwd.Contains("_rb5")){
    nrbx = 5;
  }
  else if(pwd.Contains("_rb7")){
    nrbx = 7;
  }

  printf("check nrbx rebin %d\n", nrbx);

  //original.root
  TFile *fhistroot = new TFile(argv[1]);
  if(!fhistroot->IsOpen()){
    printf("no root!!!!\n");
    exit(1);
  }

  fmcmom =(TH2D*) fhistroot->Get("mcmom");
  if(!fmcmom){
    printf("drawpar no fmcmom!\n");gDirectory->ls(); exit(1);
  }

  fmcpt =(TH2D*) fhistroot->Get("mcpt");
  if(!fmcpt){
    printf("drawpar no fmcpt!\n");gDirectory->ls(); exit(1);
  }
  
  fmcz =(TH2D*) fhistroot->Get("mcz");
  if(!fmcz){
    printf("drawpar no fmcz!\n");gDirectory->ls(); exit(1);
  }
  
  fmcxi =(TH2D*) fhistroot->Get("mcxi");
  if(!fmcxi){
    printf("drawpar no fmcxi!\n");gDirectory->ls(); exit(1);
  } 

  TH2D * hh[]={fmcz, fmcxi};
  const Int_t hn = sizeof(hh)/sizeof(TH2D*);
  for(Int_t ima=0; ima<hn; ima++){
    //rebin and then divide
    for(Int_t iy=1; iy<= hh[ima]->GetNbinsY(); iy++){
      for(Int_t ix=1; ix<= hh[ima]->GetNbinsX(); ix+=nrbx){
        const Double_t cont = hh[ima]->Integral(ix, ix+nrbx-1, iy, iy)/nrbx;
        for(Int_t ib=0; ib<nrbx; ib++){
          hh[ima]->SetBinContent(ix+ib,iy, cont);
        }
      }
    }
  }


  getMCETAall();

  printf("getMCETAalldone!\n");
  return 0;
}
